Bounding the mass of the graviton with 
gravitational waves: Effect of higher harmonics in 
gravitational waveform templates 

K. G. and Clifford M. Will^'^ 

E-mail: arun@physics.wustl.edu, cmw@wuphys.wustl.edu 

^ GRcCO, Institut d'Astrophysiquc dc Paris, UMR 7095-CNRS, Universite Pierre et 

Marie Curie, 98''*'* Bd. Arago, 75014 Paris, France 

2 LAL, Universite Paris Sud, IN2P3/CNRS, Orsay, France 

McDonnell Center for the Space Sciences, Department of Physics, Washington 
University, St. Louis MO 63130 USA 



Abstract. Observations by laser interferometric detectors of gravitational waves 
from inspiraling compact binary systems can be used to search for a dependence of the 
waves' propagation speed on wavelength, and thereby to bound the mass or Compton 
wavelength of a putative graviton. We study the effect of including higher harmonics, 
as well as their post-Newtonian amplitude corrections, in the template gravitational 
waveforms employed in the process of parameter estimation using matched filtering. 
We consider the bounds that could be achieved using advanced LIGO, a proposed third 
generation instrument called Einstein Telescope, and the proposed space interferometer 
LISA. We find that in all cases, the bounds on the graviton Compton wavelength 
are improved by almost an order of magnitude for higher masses when amplitude 
corrections are included. 
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1. Introduction and summary 

The ongoing development of laser interferometric gravitational-wave observatories on 
the ground and in space brings closer the day when gravitational radiation will be used 
as a tool for astronomical discovery and for testing fundamental physics [1] . The Laser 
Interferometer Gravitational- Wave Observatory (LIGO) [2] recently finished its fifth 
science run, operating for a year and half at its initial design sensitivity. Its European 
counterpart Virgo [3] ran for four months in coincidence with the final period of LlGO's 
science run. Both detectors will undergo a series of upgrades to improve their sensitivity 
by a factor of about ten and to be back on the air in the 2014 time frame. A third 
generation European interferometer provisionally called Einstein Telescope (ET) |1] is 
being planned to have unprecedented low-frequency sensitivity, with a seismic cutoff 
close to 1 Hz. While the ground-based detectors will be sensitive to high-frequency 
gravitational waves (1 — 10^ Hz), the proposed Laser Interferometer Space Antenna 
(LISA) [5], will be sensitive to low- frequency gravitational waves in the milli-Hertz range. 
LISA will therefore be able to detect gravitational waves from sources like supermassive 
black hole (SMBH) binaries and will complement the ground-based detectors which are 
to be sensitive to stellar mass and intermediate mass black hole (BH) binaries. Since 
the gravitational waves from these compact binaries can be very precisely modelled 
using analytical and numerical relativity, the detection and parameter estimation will 
be performed by the technique of matched filtering, whereby theoretically generated 
waveforms are used as templates to search for gravitational-wave signals in the data. 

In previous papers [6], [3, [8], [9] we studied the extent to which observations of such 
gravitational waves could test gravitational theory, in particular test whether the speed 
of gravitational waves depends on their wavelength, or, to use a shorthand phrase, 
whether the graviton has a mass (the waves themselves are completely describable by 
non-quantum physics, of course). We found that the upgraded ground-based detectors, 
provisionally denoted by AdvLIGO and AdvVirgo, could place a lower bound on the 
graviton Compton wavelength comparable to bounds (~ 10^^ km) obtained from solar 
system dynamics, while LISA could do some four orders of magnitude better. 

The basic idea is simple: if there is a mass associated with the propagation of 
gravitational waves ("a massive graviton"), then the speed of propagation will depend 
on wavelength in the form Vg ^ 1 — (A/Ag)^, where Xg is the Compton wavelength of the 
graviton, in the limit where A ^ Ag. Irrespective of the nature of the alternative theory 
that predicts a massive graviton (and notwithstanding the difficulties in defining such 
theories free of pathologies such as the ZvDV discontinuity [lOl [H]), it is reasonable to 
expect the differences between such a hypothetical theory and general relativity in the 
predictions for the evolution of massive compact binaries to be of order (A/A^)^, and 
therefore to be very small, given that A ~ 10^ km for stellar mass inspirals and ~ 10^ 
km for massive black hole inspirals. 

As a result, the gravitational waveform seen by an observer close to the source will 
be very close to that predicted by general relativity. However, as seen by a detector 
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at a distance D, hundreds to thousands of Mpc away, the phasing of the signal will be 
distorted because of the shifted times of arrival, At ~ D[\/\gY of waves emitted with 
different wavelengths during the inspiral. In addition to measuring the astrophysical 
parameters of the system, such as masses and spins, the matched filtering technique 
permits one to estimate or bound such effects. We point out that small deviations from 
general relativity will not strongly affect the detection of the gravitational wave signals. 
The reason is that, for detection, one maximizes the overlap function of the signal with 
templates over the parameters of the template (see [I2] for details). Even if the actual 
signal is slightly different from that predicted by GR, templates based on GR should 
be able to recover it but with significant biases in the estimated parameters. Parameter 
estimation hence should be seen as a post-detection problem. 

In our earlier matched filtering analyses [HI [7l [HI E], we chose a particular form 
for the theoretical waveform template, known as the "restricted waveform (RWF)", 
constructed using the dominant quadrupole amplitude of the wave evaluated in the 
lowest-order "Newtonian" approximation, but with a phase expressed to the highest 
post-Newtonian (PN) order available at the time. This was generally second post- 
Newtonian (2PN) order, or 0{\v/cf) beyond the leading quadruple approximation, and 
included the effects of non-precessing spin (see [13] for a review of the post-Newtonian 
phasing formulae). The Fourier domain waveform in this approximation has only the 
leading order amplitude multiplying a term whose phase is proportional to twice the 
orbital phase. 

However, recent work has pointed out that the inclusion of other multipoles in the 
wave amplitude, as well as PN corrections to the amplitudes, can have a dramatic effect 
on the estimation of certain system parameters [MfT5lfM[T7l[T8l[T9|[2ni[2T|[22l[23| 
[2H125]. The effects of incorporating higher harmonics include increased mass reach for 
the ground-based [19] and space-based [21] detectors, improved estimation of the mass 
parameters [201 [22l [231 El] , improved angular resolution [161 ESI ESI El] and improved 
estimation of cosmological parameters such as those associated with the dark energy 
equation of state [221 125] . 

These improvements arise from two effects, (i) The presence of higher harmonics 
of the orbital phase increases the frequency span of the signal in the detector band, (ii) 
The structure of the waveform is richer because, even though they are smaller than the 
dominant quadrupole term, the new amplitude terms and their PN corrections introduce 
new functions of masses, spins, and inclination angles. 

The purpose of this paper is to explore the effects of using the full waveform (FWF) 
on the bounds on A^. In particular, we revisit the earlier bounds for non-spinning 
binaries by including the effects of PN amplitudes up to 3PN order [26l [271 [23 [29] and 
phasing up to 3.5 PN order [301 [SD [321 [23], including harmonics as high as 8 times the 
orbital phase. We also consider the three different detector configurations mentioned 
above: the second generation AdvLIGO, the proposed third generation ET, and the 
proposed space-based interferometer LISA. 

Fig. [1] displays our central results. The l-cr bounds on \g obtainable from AdvLIGO, 
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Figure 1. Bounds on the graviton Compton wavelength that can be deduced from 
AdvLIGO, Einstein Telescope and LISA. The mass ratio is 2. The distance to the 
source is assumed to be 100 Mpc for AdvLIGO and ET, and 3 Gpc for LISA. 



ET and LISA are plotted as a function of the total mass of the binary for a fixed mass 
ratio of m2/mi = 2. For AdvLIGO and ET, the source is assumed to be at a luminosity 
distance of 100 Mpc and for LISA the SMBH binary is assumed to be 3 Gpc away. 
The bounds from the Newtonian RWF and 3PN FWF are compared. Inclusion of 
amplitude corrections and the higher harmonics improve the bounds for both ground- 
based configurations and at the high-mass end for LISA. The improvement is more 
than an order of magnitude for heavier binaries, because higher harmonics play a more 
prominent role for such systems. Typical bounds, with the use of higher harmonics, 
for AdvLIGO, ET and LISA are 10^^ km, 10^^ km and lO'^^ km, respectively. The best 
bound, not surprisingly, will be provided by LISA, thanks to its low frequency sensitivity, 
to the high signal-to-noise ratios with which it will be observing the supermassive binary 
black hole coalescences, and to the very large distances involved. Though our results 
are for a specific location and orientation of the binary, we have verified that the bounds 
are not significantly altered by different source positions and orientations. 



The remainder of the paper provides details underlying these results. In Sec. [21 we 
describe the full-waveform model used, the noise curves for the various detectors, and 
the technique of matched filtering. Section [3] details the bounds obtainable from the 
various detectors. 
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2. Parameter estimation using full waveform templates 

As our waveform model we begin with amplitude-corrected, general relativistic 
waveforms which are 3PN accurate in amplitude and 3.5PN accurate in phasing. We 
ignore the spins of the bodies in the binary system. Previous calculations used waveforms 
which are of Newtonian order in amplitude and 2PN order in phase. As opposed to the 
Newtonian waveforms, the 3PN amplitude-corrected waveforms contain all harmonics 
from \& up to 8 where ^ is the orbital phase (the leading quadrupole component is 
at 2*). 

The effect of a massive graviton is included in the expression for the orbital phase 
following Ref. [6]. The wavelength-dependent propagation speed changes the arrival 
time ta of a wave of a given emitted frequency /e relative to that for a signal that 
propagates at the speed of light; that time is given, modulo constants,by 

—1 (1) 

where fe and are the wave frequency and time of emission as measured at the emitter, 
respectively, Z is the cosmo logical redshift, and 

n = il±^ a{t)dt , (2) 

'^0 Jte 

where Oq = a{ta) is the present value of the scale factor (note that D is not exactly the 
luminosity distance 0). This affects the phase of the wave accordingly. In the frequency 
domain, this adds a term to the phase ip{f) of the Fourier transform of the waveform 
given by Aip^f) = —TTD/feX^. Then, for each harmonic of the waveform with index k, 
one adds the term 

AMf) = |a^(2//A:) = -^irD/fAl ■ (3) 

Here k = 2 denotes the dominant quadrupole term, with phase 2\l/, k = 1 denotes the 
term with phase = 3 denotes the term with phase 3\Ef, and so on. 

This is an adhoc procedure because a massive graviton theory will undoubtedly 
deviate from GR not just in the propagation effect, but also in the way gravitational wave 
damping affects the phase, as well as in in the amplitudes of the gravitational waveform. 
If, for example, such a theory introduces a leading correction to the quadrupole phasing 
'ipqua.d ~ (ttTW/e)"^''^ of Order (A/Ag)^x (vrTW/e)"^/^, where M is the chirp mass, then the 
propagation induced phasing term ([3]) will be larger than this correction term by a factor 
of order k'^{D / M.){'KM.feT^^ ~ {D/M)v^. Since v ~ 0.1 for the important part of the 
binary inspiral, and D ~ hundreds to thousands of Mpc, it is clear that the propagation 
term will dominate. In any case, given the fact that there is no generic theory of a 
massive graviton, we have no choice but to omit these unknown contributions. 

I For Z ^ 1, D is roughly equal to luminosity distance D^. Hence we have assumed D ~ Dl in the 
case of ground based detectors for which we consider sources at 100 Mpc. For LISA, we have carefully 
accounted for this difference. 
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However the effect of our assumption of neglecting tfie spin of the binary could 
be more severe. As studied in detail in [8], the spin effects could weaken the bounds 
on Xg. This conclusion was based on the restricted waveforms and on nonprecessing 
spins. When the higher harmonics from the polarization as well as spin precession 
are included, the trends may be different. But this will require a thorough analysis 
of parameter estimation with precessing binary waveforms; this will be the subject of 
future work. 

2.1. Instrument noise models 

The matched filtering procedure for estimating errors requires knowledge of the noise 
characteristics of the detectors. In the frequency domain, these characteristics are 
embodied in the noise power spectral density (PSD). We use analytical fits to the 
designed PSDs of various detector configurations. For AdvLlGO, we have used the 
analytical fit provided in Eq. (3.7) of [M]. We assume the ET to be an L-shaped detector 
with an arm length of 10 km. The PSD we use is from Ref. [35]. This is essentially 
the analytical fit given in Eqs (4.4) and (4.5) of Ref. [20] with slight modifications to 
incorporate the arm length. 

For the case of LISA, we have used the noise PSD from Ref. |8], Eqs (2.28)-(2.32), 
though there have been minor modifications to the LISA noise PSD since then. This 
choice of the noise curve helps us to calibrate our codes by reproducing the results of 
Ref. [8] for the restricted waveform case. 

2.2. Frequency cut-offs 

The gravitational wave signals are converted from the time domain to the frequency 
domain using the stationary-phase approximation (SPA). We assume the signals to last 
in the respective detector bands between frequencies /lower and /upper- We assume that 
the upper cut-off is given by the last stable orbit (LSO) of the system. Hence the /c— th 
harmonic, for instance, is sharply cut-off at k -Flso using step functions, where -Flso is 
the orbital frequency at the last stable orbit. Readers can refer to Refs. [19l EH [36] 
for the details and subtleties associated with this method and related assumptions. See 
Refs [371 IM] for a general discussion of the step function cut-offs employed in SPA. 

For the ground-based detectors, the low frequency cut-off is fixed by the seismic 
cut-off of the detector. Following the designed noise PSDs, we assume this to be 20 Hz 
for AdvLIGO and 1 Hz and 10 Hz for Einstein Telescope. 

For the case of LISA, the low frequency cut-off is defined in a different way. Since 
LISA can observe the sources for extended periods of order months to years, we make 
the standard assumption that the source is observed in the LISA band for one year 
prior to merger. For a particular SMBH system, we use -Flso to be a measure of the 
upper cut-off orbital frequency and deduce the low orbital frequency cut-off assuming a 
1 year (or less) duration of the signal in the LISA band. Hence both the low frequency 
and high frequency cut-offs will be different for each harmonic (see Ref. [21] for details). 
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Following Ref. [5], we choose 

/lower = Max (10"^HZ, /be, 
/begin = 4.149 X 10-^ f-J 



M 



'gin 



) , where 
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3/8 



(5) 



(4) 



where Ai is the chirp mass of the binary and Tobs is the duration of observation of the 
signal, which is assumed to be 1 year (or less). Notice also that we have assumed, rather 
conservatively, that LISA is not sensitive to frequencies below 10"^ Hz. In our choice 
the k-th harmonic will last from fc/iowcr/2 until fc-FLSo- 

2.3. Assumptions about the sources 

We have assumed the typical distance to the binary black holes to be 100 Mpc for the 
AdvLIGO and ET cases, and 3 Gpc for LISA. Notice, however, that the bounds on 
Xg are roughly independent of the distance of the source for a given detector [6]. This 
is because the size of the massive graviton term in the phasing, Eq. ([3]), is directly 
proportional to the distance of the source, while the errors in estimating or bounding 
this term are roughly inversely proportional to signal-to-noise ratio, and therefore also 
increase linearly with distance. 

All the plots which we provide assume a particular source position and orientation 
of the binary. These fixed source coordinates are specified with respect to an earth-based 
coordinate system for the ground based detectors, while for LISA they are specified with 
respect to the fixed LISA frame (we ignore LISA's orbital motion). 

The bounds we derive will vary with the source's location and orientation for all 
three detector configurations. We performed a numerical experiment to estimate crudely 
this uncertainty in our estimate by deriving estimates for 100 random realizations of the 
source direction. We found, for example, that the bound on varied between 1 and 
8 X 10^^ km for AdvLIGO for a AOMq — 80Mq system. Hence we expect the typical 
orders of magnitude we quote to be representative of the bounds obtainable on average. 
Note that all the bounds that we quote are for single GW events. In the happy event 
of numerous detections, then the bounds could improve either in value or in confidence 
level; however we have not analyzed this possibility. 

2.4. Errors obtained using the Fisher matrix approach 

Our estimate of the bounds on the massive graviton parameter is based on the Fisher 
matrix formalism. We construct the Fisher matrix for the different detector noise PSDs 
using the amplitude corrected PN waveform model described earlier, converted to the 
Fourier domain using the stationary phase approximation. We use a six- dimensional 
parameter space consisting of the time and phase {tc, (j)c) of coalescence, the chirp mass 
A4, the mass asymmetry parameter S = |mi — m2|/(mi -|- m2), the massive graviton 
parameter Pg = 7r^D7W/A^(l + Z), and the luminosity distance D^. We fix the three 
angles, 6, and ip which appear in the antenna pattern functions to be tt/S, tt/G and 
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Figure 2. Bounds on Ag using restricted and 3PN full waveforms for AdvLIGO (left 
panel) and Einstein Teleseope (right panel). For Einstein Telescope, we show the 
bounds assuming two possible seismic cut-off frequencies: 10 Hz and 1 Hz. 

7r/4 respectively and the inclination angle of the binary to be i = 7r/3. Details of the 
Fisher matrix approach as applied to the compact binary coalescence signals can be 
found in Refs. [391 SOI El], and more recently in Ref. [H] ,which critically reexamines 
the caveats involved in using the Fisher matrix formalism to deduce error bounds for 
various gravitational wave detector configurations. 

The square root of each of the diagonal entries in the inverse of the Fisher matrix 
gives a lower bound on the error covariance of any unbiased estimator. Our focus here 
is solely on the diagonal element corresponding to the massive graviton parameter. The 
1 — cr error bar on fig can be translated into a bound on the Compton wavelength using 
A/3g = (3g, and this is the quantity that we use in the plots as well as in the discussions. 

In the next section we discuss our results in detail for the three detector 
configurations considered. 

3. Bounds on massive graviton theories from AdvLIGO, ET and LISA 
3.1. Advanced LIGO and Einstein Telescope 

Figure [2] shows the bounds on Xg possible from future GW observations with 
AdvLIGO and Einstein Telescope. There is evident improvement from using the FWF 
in the AdvLIGO case. For binaries which are well-detected with the RWF, the inclusion 
of higher harmonics improves the massive graviton bounds by a factor of a few. The 
most striking feature is that even for masses that are beyond the range detectable by 
RWF templates, the FWF based templates still put bounds which are better than the 
best bounds from RWF. The best bounds would be for intermediate mass BH binaries 
whose total mass lies in the range 50-100 Mq. 

Because of its improved sensitivity, ET will be able to put more stringent bounds, 
roughly an order of magnitude better, compared to AdvLIGO. On the right panel 
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Figure 3. Bound on Xg using LISA, as a function of total mass for restricted and full 
waveforms. Results for different mass ratios 5 -100 are shown. 



of Fig. [21 we show the bounds that are possible with two different seismic cut-off 
frequencies, 10 Hz and 1 Hz, for ET. The FWF results for a 10 Hz cut-off are as good as 
the results for RWF with a 1 Hz cut-off. This is because including the higher harmonics 
produces roughly the same effect as having a lower cut-off frequency because both extend 
the sensitivity of the detector to lower frequencies. However the best estimates for ET 
are still from the case where the seismic cut-off is 1 Hz and FWF templates are employed 
and for intermediate mass BH binaries whose masses are of the order of a few hundreds 
of solar masses. Also it is interesting to note that templates which account for higher 
harmonics essentially bridge the gap in mass coverage of different gravitational wave 
detectors (see Fig. [T]). The "steps" in the curves for AdvLIGO and ET at the high mass 
ends are due to the effects on different harmonics of the sharp step function cut-offs at 
high frequencies. 

In the FWF cases for both AdvLIGO and ET we have checked, for all mass ranges 
shown in the plot, that the Fisher matrix is not ill-conditioned and the binaries are 
observed with an SNR of at least 8. There were cases for which the SNR was less than 
5 but the Fisher matrix was found to be well-conditioned, but we have not used those 
systems because for very low SNRs, the Fisher matrix formalism itself breaks down |41j . 

3.2. LISA 

In Fig. 3 we show the possible bounds on from observations by LISA of supermassive 
binary inspirals for different mass ratios. The bounds are almost three orders 
of magnitude better than those from the best ground-based experiments. This 
improvement is almost entirely due to the larger masses involved, as can be seen from 
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the dependence of the bound on Xg on the relevant parameters of the system and the 
detectors, given by Eq. (4.9) of [6]: 



characteristic "knee" frequency, or frequency where the PSD is a minimum. The 
quantities /(7) and A are determined from the Fisher matrix inversion and are largely 
independent of either 5*0 or /q, or of the SNR of the signal. In any case, the bound is 
only weakly dependent on these variables. The ratio D/{1 + Z)Dl is weakly dependent 
on distance, reflecting the fact that the effect of the massive graviton and the estimation 
errors both grow with distance. As it turns out, the factor S^^ fl^'^ is roughly the same 
for LISA as it is for AdvLIGO, and thus the bound on is almost entirely proportional 
to the chirp mass of the source. 

Though the bound decreases with increase in the mass ratio, the improvement due 
to the inclusion of higher harmonics is more profound with increasing mass ratio. The 
amplitude corrections to the higher harmonics are proportional to the mass asymmetry 
of the binary for all odd PN amplitude orders (see e.g. Eqs (8.9)-(8.10) of [29]) which 
include the leading correction at 0.5PN order, and hence the above observation is not 
surprising. It should also be borne in mind that, for very large mass ratios, the systems 
are more like Extreme Mass Ratio Inspirals, and many other PN effects will have to be 
accounted for in order to detect them. As in the case of the ground-based detectors, 
use of the FWF increases the range of masses over which stringent bounds can be put 
on the Xg term. This range of masses (IO^Mq- a few times IO'^'Mq) is very much in the 
astrophysically interesting regime. 

As mentioned earlier, we have not taken into account the orbital motion of LISA 
in deriving these estimates. However, the addition of angular parameters should not 
strongly affect the reported bounds because they are rather uncorrelated with the (3g 
parameter. On the contrary, inclusion of LISA's orbital motion could enhance the SNRs 
and thereby improve the bounds. 

4. Conclusions 

We have shown that the use of amplitude corrected PN waveforms which incorporate 
higher harmonics of the orbital phase in constructing templates for matched filtering 
(especially for parameter estimation) can provide more stringent bounds on the massive 
graviton parameter. Third generation ground-based interferometers such an Einstein 
telescope, and the space-based LISA would almost continuously cover the BH inspirals 
ranging from stellar mass binaries of a few tens of solar mass up to supermassive BHs 



Table 1 summarizes the typical bounds on the graviton Compton wavelength 
achievable from different interferometric gravitational-wave detectors incorporating 





where 5*0 is a parameter that establishes the scale of the PSD (in Hz 



), /o is a 



of WM, 
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Detector 


Mass Range (M©) 


SNR 


Xg (km) 


AdvLIGO 


1-400 


few tens 


10^2 


Einstein Telescope 


1-103-4 


few tens-few hundreds 


1013-14 


LISA 




few thousands 


10^6 



Table 1. Summary of bounds on the graviton Compton wavelength from various 
gravitational wave detectors using full waveform templates 



higher harmonics into the gravitational waveform templates. Also shown are typical 
orders of magnitudes of the observational mass range and the signal to noise ratio 
(SNR). The mass range varies with the chosen low frequency cut-off of the detector. 
For ET, the mass range is between IMq and IO^Mq or IO'^Mq, depending on whether 
the seismic cut-off chosen is 10 Hz or 1 Hz, respectively. For LISA, the mass range is 
between IO^Mq and IO^Mq or 10^ Mq, depending on whether the low frequency cut-off 
of 10"^ or 10"^ Hz, respectively . Numbers for the mass range and SNR have been taken 
from the literature [191 El]- The best bounds would be from the proposed space-based 
interferometer, LISA. 
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